# Set working directory to main replication folder
setwd("C:/Users/Eric/Dropbox/Book/analysis/replication/")

# Load libraries
library(plyr)
library(dplyr)
library(tidyverse)
library(ggplot2)
library(lubridate)
library(MASS)
library(zoo)
library(AER)
library(broom)
library(data.table)
library(extrafont)
library(foreach)
library(doParallel)
library(lmtest)
library(scales)
library(stargazer)

# Load functions
source("./function_code/battle_coef_plot.R")
source("./function_code/msm_coef_plot.R")
source("./function_code/msm_coef_table.R")
source("./function_code/combine_rev_results.R")

# Read in data
d = fread("./data/wardata_book.csv")
d = d |> arrange(warName, dayNum)

# Negotiation start data
s = fread("./data/negstarts.csv")


#### A QUIETER BATTLEFIELD  ####
controlsB = "ceasefire + AOS60 + AOA + issues + contig + cincP + polBinA + maj + nuk + nBellig + keyAllies + anyDipEnemy + logBattles + nBattleTrend"

batFit1a = glm(as.formula(paste("nBattle ~ neg + factor(warName)",
                                controlsB, sep="+")), data=d, family="poisson")
batFit1ar = coeftest(batFit1a, vcov=vcovCL, cluster = ~ warName)

batFit1b = glm(as.formula(paste("nBattle ~ negInt + negExt + factor(warName)",
                                controlsB, sep="+")), data=d, family = "poisson")
batFit1br = coeftest(batFit1b, vcov=vcovCL, cluster = ~ warName)


# Change in number of battles
round(exp(batFit1a$coefficients["neg"]), 3)    # 0.699
round(exp(batFit1b$coefficients["negExt"]), 3) # 0.943
round(exp(batFit1b$coefficients["negInt"]), 3) # 0.57

car::linearHypothesis(batFit1b, "negInt = negExt", singular.ok=T) # p << 0.01


# Overdispersion?
round(mean(d$nBattle), 3) # 1.038
round(var(d$nBattle), 3)  # 6.142
dispersiontest(batFit1a)  # p = 1
dispersiontest(batFit1b)  # p = 1


#### FIGURE 4.1: Coefficient plots for Poisson models of active battles with war fixed effects (n = 36,834) ####
battleCoefPlot(batFit1ar, batFit1br, "Main")
file.rename("./figures/battleCoefPlotMain.pdf",
            "./figures/figure_4.1.pdf")


# Proportion of negotiations with ceasefires 
sum(d$neg==1 & d$ceasefire==1)/sum(d$neg==1) # 0.091

# Proportion of ceasefires with negotiations
sum(d$neg==1 & d$ceasefire==1)/sum(d$ceasefire==1) # 0.300

#### Measure similarity of negotiation and ceasefire measures ####
jaccard <- function(a, b) {
  intersection = length(intersect(a, b))
  union = length(a) + length(b) - intersection
  return (intersection/union)
}

jaccard(d$neg, d$ceasefire) # Essentially 0


#### CHANGING BATTLEFIELD TRENDS ####

s0 = s |> filter(endwar != "1")  # All failed negotiations
s1 = s |> filter(endwar != "1") |> filter(resNum == 1)  # Only first failed negotiations per war 
s2 = s |> filter(endwar != "1") |> filter(resNum > 1)   # Only subsequent failed negotiations per war

slist = list(s0, s1, s2)
slabels = c("all", "early", "later")
fights = c("AOS60", "outSum60", "apropdiff", "apropdiffDOE")
bands = c(30, 60, 90)

revList = tidyRevList = robustList = list()

for (q in 1:length(slist)) {
  ss = slist[[q]]
  
  resList = tidyResList = stargazerList = robList = list()
  
  for (f in 1:length(fights)) {
    
    fight = fights[f]
  
    
    coefLista = coefListb = tidyLista = tidyListb = robLista = robListb = list()
    
    for (b in 1:length(bands)) {
      band = bands[b]
      
      negList = list()
      for (i in 1:nrow(ss)) {
        oneneg = ss[i,]
        nstart = as.Date(oneneg$startdate)
        nend = as.Date(oneneg$enddate)
        
        # Get the war info
        onewar = d |> filter(warName==oneneg$warName)
        
        # Get the days before 
        befstart = max(nstart - band, min(onewar$date))
        befDays = onewar[which(onewar$date %in% as.Date(befstart:(nstart-1))),]
        befDays$afterNeg = 0
        
        # Get the days after
        aftend = min(nend + band, max(onewar$date))
        afterDays = onewar[which(onewar$date %in% as.Date((nend+1):aftend)),]
        afterDays$afterNeg = 1
        
        baneg = rbind(befDays, afterDays)
        
        negList[[i]] = baneg
      }
      pp3 = do.call("rbind", negList)
      pp3 = as.data.frame(pp3)
      
      resPostFit1a = lm(pp3[,fight] ~ afterNeg + factor(warName), data=pp3)
      resPostFit1at = tidy(resPostFit1a) |> filter(!grepl("factor", term))
      coefLista[[b]] = resPostFit1a
      tidyLista[[b]] = resPostFit1at
      robLista[[b]] = coeftest(resPostFit1a, vcov=vcovHC, cluster = ~ warName)
      
      resPostFit1b = lm(pp3[,fight] ~ afterNeg + neg + dayNum + factor(warName), data=pp3)
      resPostFit1bt = tidy(resPostFit1b) |> filter(!grepl("factor", term))
      coefListb[[b]] = resPostFit1b
      tidyListb[[b]] = resPostFit1bt
      robListb[[b]] = coeftest(resPostFit1b, vcov=vcovHC, cluster = ~ warName)
    }
    
    names(coefListb) = paste0("days", bands)
    names(tidyListb) = paste0("days", bands)
    names(robListb) = paste0("days", bands)
    
    # Store the statistical results
    resList[[f]] = coefListb
    tidyResList[[f]] = tidyListb
    robList[[f]] = robListb
  }
  
  names(resList) = fights
  names(tidyResList) = fights
  names(robList) = fights
  
  revList[[q]] = resList
  tidyRevList[[q]] = tidyResList
  robustList[[q]] = robList
  
}
names(revList) = slabels
names(tidyRevList) = slabels
names(robustList) = slabels


# Gather results
turn_all_aos_60 = revList$all$AOS60$days60
turn_all_mom_60 = revList$all$outSum60$days60
turn_all_apd_60 = revList$all$apropdiff$days60
turn_all_apddoe_60 = revList$all$apropdiffDOE$days60

turn_early_aos_60 = revList$early$AOS60$days60
turn_early_mom_60 = revList$early$outSum60$days60
turn_early_apd_60 = revList$early$apropdiff$days60
turn_early_apddoe_60 = revList$early$apropdiffDOE$days60

turn_later_aos_60 = revList$later$AOS60$days60
turn_later_mom_60 = revList$later$outSum60$days60
turn_later_apd_60 = revList$later$apropdiff$days60
turn_later_apddoe_60 = revList$later$apropdiffDOE$days60


# The estimated coefficient for effect of post-negotiation on recent imbalance
round(coef(turn_all_aos_60)['afterNeg'], 3) # 0.225

# The estimated coefficient for effect of post-negotiation on recent momentum
round(coef(turn_all_mom_60)['afterNeg'], 3) # -0.232

# The estimated coefficient for effect of post-negotiation on inconsistency
round(coef(turn_all_apd_60)['afterNeg'], 3) # -0.013

# Variance in these measures within wars
fightVars = d |> dplyr::group_by(warName) |>
  dplyr::summarize(varsAOS = var(AOS60, na.rm=T),
                   varsOutSum = var(outSum60, na.rm=T),
                   varsAPD = var(apropdiff, na.rm=T),
                   varsAPDDOE = var(apropdiffDOE, na.rm=T))

round(median(fightVars$varsAOS), 3)             # 0.594
round(median(fightVars$varsOutSum), 3)          # 0.805
round(median(fightVars$varsAPD, na.rm=T), 3)    # 0.005
round(median(fightVars$varsAPDDOE, na.rm=T), 3) # 0.004

# Median value of inconsistency on the last day of war
d |> dplyr::filter(terminate==1) |> summarize(mean(apropdiff, na.rm=T)) # 0.318


# Combine results for plot
rev_aos_60 = combineResults(fightMeasure = "aos", window=60)
rev_mom_60 = combineResults(fightMeasure = "mom", window=60)
rev_apd_60 = combineResults(fightMeasure = "apd", window=60) 
rev_apddoe_60 = combineResults(fightMeasure = "apddoe", window=60) 
revData = bind_rows(rev_aos_60, rev_mom_60, rev_apd_60, rev_apddoe_60)

revData$vars = factor(revData$vars, levels=c("afterNeg", "neg", "dayNum"))
revData$model = factor(revData$model, levels=c("all", "early", "later"), labels=c("All", "First", "Subsequent"))
revData$fight = factor(revData$fight, levels=c("aos", "mom", "apd", "apddoe"), labels=c("Recent imbalance", "Momentum", "Inconsistency", "Inconsistency (DOE)"))
revData$signif = factor(revData$signif, levels=c("No", "Yes"))


#### FIGURE 4.2: Coefficient plots for OLS models of changes in battlefield outcomes after failed negotiations, using a sixty-day temporal window ####
revData |> filter(fight!="Inconsistency (DOE)") |>
  ggplot(aes(est, fct_rev(vars), group=fct_rev(model))) + 
  geom_vline(xintercept=0, linetype=2) + 
  geom_pointrange(aes(xmin=lower95, xmax=upper95, color=model, shape=signif), position=position_dodge(width=0.6)) + 
  theme_bw() + xlab("Coefficient estimate") +
  scale_y_discrete(name="", labels=c("Time trend", "Concurrent\nnegotiation", "Post-\nnegotiation")) + 
  scale_x_continuous(breaks=breaks_pretty()) + 
  scale_color_manual("Negotiations", 
                     values=c("black", "gray30", "gray60")) +
  scale_shape_manual("Statistically\nsignificant (95%)", values=c(16, 15)) +
  facet_grid(cols=vars(fight), scales="free") + 
  theme(panel.spacing.x = unit(6, "mm"),
        legend.position = "bottom",
        legend.direction = "vertical")

ggsave("./figures/figure_4.2.pdf", width=7, height=4.5)


######################################################################################################


#### APPENDIX MATERIALS FOR CHAPTER 4 ####

##### APPENDIX 4.1: FULL POISSON REGRESSION RESULTS #####

###### TABLE A16: Results of Poisson regressions ######

# Initial results before clustered standard errors
# Note that the number of observations per model only shows up here
stargazer(batFit1a, batFit1b,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum", "AOS60:logDayNum", "warName", "issues", "contig"),
          order = c("neg", "negInt", "negExt", 
                    "ceasefire", "AOS60", "AOA", "issues", "contig", "cincP", "polBinA",
                    "maj", "nuk", "nBellig", "anyDipEnemy", "keyAllies", "logBattles"),
          covariate.labels = c("All negotiation", "Internal negotiation", "External negotiation",
                               "Ceasefire", "Recent imbalance", "Overall imbalance",
                               "CINC ratio",
                               "Democratic initiator", "Major power", "Nuclear", "Number of states",
                               "Opp. dip. representation",
                               "Major allies",
                               "Completed battles", "Active battle trend"),
          dep.var.labels.include=T,
          dep.var.labels = c("Active battles", "Active battles")
)

# Results after clustered standard errors
# Note that the number of observations per model does not show up here
stargazer(batFit1ar, batFit1br,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum", "AOS60:logDayNum", "warName", "issues", "contig"),
          order = c("neg", "negInt", "negExt", 
                    "ceasefire", "AOS60", "AOA", "issues", "contig", "cincP", "polBinA",
                    "maj", "nuk", "nBellig", "anyDipEnemy", "keyAllies", "logBattles"),
          covariate.labels = c("All negotiation", "Internal negotiation", "External negotiation",
                               "Ceasefire", "Recent imbalance", "Overall imbalance",
                               "CINC ratio",
                               "Democratic initiator", "Major power", "Nuclear", "Number of states",
                               "Opp. dip. representation",
                               "Major allies",
                               "Completed battles", "Active battle trend"),
          dep.var.labels.include=T,
          dep.var.labels = c("Active battles", "Active battles")
)




##### APPENDIX 4.2: CHANGING BATTLEFIELD TRENDS #####

###### APPENDIX 4.2.1: Full 60-Day Results ######

####### TABLE A17: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using all negotiations and a 60-day temporal window around failed negotiations 
stargazer(turn_all_aos_60, turn_all_mom_60, turn_all_apd_60,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

####### TABLE A18: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using the first negotiation in each war and a 60-day temporal window around failed negotiations 
stargazer(turn_early_aos_60, turn_early_mom_60, turn_early_apd_60,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

####### TABLE A19: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using all negotiations after the first in each war and a 60-day temporal window around failed negotiations
stargazer(turn_later_aos_60, turn_later_mom_60, turn_later_apd_60,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))


###### APPENDIX 4.2.2: Full 30-Day Results ######

# Gather results
turn_all_aos_30 = revList$all$AOS60$days30
turn_all_mom_30 = revList$all$outSum60$days30
turn_all_apd_30 = revList$all$apropdiff$days30
turn_all_apddoe_30 = revList$all$apropdiffDOE$days30

turn_early_aos_30 = revList$early$AOS60$days30
turn_early_mom_30 = revList$early$outSum60$days30
turn_early_apd_30 = revList$early$apropdiff$days30
turn_early_apddoe_30 = revList$early$apropdiffDOE$days30

turn_later_aos_30 = revList$later$AOS60$days30
turn_later_mom_30 = revList$later$outSum60$days30
turn_later_apd_30 = revList$later$apropdiff$days30
turn_later_apddoe_30 = revList$later$apropdiffDOE$days30


####### TABLE A20: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using all negotiations and a 30-day temporal window around failed negotiations 
stargazer(turn_all_aos_30, turn_all_mom_30, turn_all_apd_30,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

####### TABLE A21: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using the first negotiation in each war and a 30-day temporal window around failed negotiations 
stargazer(turn_early_aos_30, turn_early_mom_30, turn_early_apd_30,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

####### TABLE A22: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using all negotiations after the first in each war and a 30-day temporal window around failed negotiations
stargazer(turn_later_aos_30, turn_later_mom_30, turn_later_apd_30,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

# Combine results for plot
rev_aos_30 = combineResults(fightMeasure = "aos", window=30)
rev_mom_30 = combineResults(fightMeasure = "mom", window=30)
rev_apd_30 = combineResults(fightMeasure = "apd", window=30) 
rev_apddoe_30 = combineResults(fightMeasure = "apddoe", window=30)
revData30 = bind_rows(rev_aos_30, rev_mom_30, rev_apd_30, rev_apddoe_30)

revData30$vars = factor(revData30$vars, levels=c("afterNeg", "neg", "dayNum"))
revData30$model = factor(revData30$model, levels=c("all", "early", "later"), labels=c("All", "First", "Subsequent"))
revData30$fight = factor(revData30$fight, levels=c("aos", "mom", "apd", "apddoe"), labels=c("Recent imbalance", "Momentum", "Inconsistency", "Inconsistency (DOE)"))
revData30$signif = factor(revData30$signif, levels=c("No", "Yes"))

####### FIGURE A10: Coefficient plots for OLS models of changes in battlefield outcomes after failed negotiations, using a 30-day temporal window #######
revData30 |> filter(fight!="Inconsistency (DOE)") |>
  ggplot(aes(est, fct_rev(vars), group=fct_rev(model))) + 
  geom_vline(xintercept=0, linetype=2) + 
  geom_pointrange(aes(xmin=lower95, xmax=upper95, color=model, shape=signif), position=position_dodge(width=0.6)) + 
  theme_bw() + xlab("Coefficient estimate") +
  scale_y_discrete(name="", labels=c("Time trend", "Concurrent\nnegotiation", "Post-\nnegotiation")) + 
  scale_x_continuous(breaks=breaks_pretty()) + 
  scale_color_manual("Negotiations", 
                     values=c("black", "gray30", "gray60")) +
  scale_shape_manual("Statistically\nsignificant (95%)", values=c(16, 15)) +
  facet_grid(cols=vars(fight), scales="free") + 
  theme(panel.spacing.x = unit(6, "mm"),
        legend.position = "bottom",
        legend.direction = "vertical")

ggsave("./figures/figure_a10.pdf", width=7, height=4.5)


###### APPENDIX 4.2.3: Full 90-Day Results ######

# Gather results
turn_all_aos_90 = revList$all$AOS60$days90
turn_all_mom_90 = revList$all$outSum60$days90
turn_all_apd_90 = revList$all$apropdiff$days90
turn_all_apddoe_90 = revList$all$apropdiffDOE$days90

turn_early_aos_90 = revList$early$AOS60$days90
turn_early_mom_90 = revList$early$outSum60$days90
turn_early_apd_90 = revList$early$apropdiff$days90
turn_early_apddoe_90 = revList$early$apropdiffDOE$days90

turn_later_aos_90 = revList$later$AOS60$days90
turn_later_mom_90 = revList$later$outSum60$days90
turn_later_apd_90 = revList$later$apropdiff$days90
turn_later_apddoe_90 = revList$later$apropdiffDOE$days90


####### TABLE A23: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using all negotiations and a 90-day temporal window around failed negotiations 
stargazer(turn_all_aos_90, turn_all_mom_90, turn_all_apd_90,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

####### TABLE A24: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using the first negotiation in each war and a 90-day temporal window around failed negotiations 
stargazer(turn_early_aos_90, turn_early_mom_90, turn_early_apd_90,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

####### TABLE A25: Results of ordinary least squares regressions analyzing changes in battlefield measures, #######
# using all negotiations after the first in each war and a 90-day temporal window around failed negotiations
stargazer(turn_later_aos_90, turn_later_mom_90, turn_later_apd_90,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))


# Combine results for plot
rev_aos_90 = combineResults(fightMeasure = "aos", window=90)
rev_mom_90 = combineResults(fightMeasure = "mom", window=90)
rev_apd_90 = combineResults(fightMeasure = "apd", window=90) 
rev_apddoe_90 = combineResults(fightMeasure = "apddoe", window=90) 
revData90 = bind_rows(rev_aos_90, rev_mom_90, rev_apd_90, rev_apddoe_90)

revData90$vars = factor(revData90$vars, levels=c("afterNeg", "neg", "dayNum"))
revData90$model = factor(revData90$model, levels=c("all", "early", "later"), labels=c("All", "First", "Subsequent"))
revData90$fight = factor(revData90$fight, levels=c("aos", "mom", "apd", "apddoe"), labels=c("Recent imbalance", "Momentum", "Inconsistency", "Inconsistency (DOE)"))
revData90$signif = factor(revData90$signif, levels=c("No", "Yes"))


####### FIGURE A11: Coefficient plots for OLS models of changes in battlefield outcomes after failed negotiations, using a 90-day temporal window #######
revData90 |> filter(fight!="Inconsistency (DOE)") |>
  ggplot(aes(est, fct_rev(vars), group=fct_rev(model))) + 
  geom_vline(xintercept=0, linetype=2) + 
  geom_pointrange(aes(xmin=lower95, xmax=upper95, color=model, shape=signif), position=position_dodge(width=0.6)) + 
  theme_bw() + xlab("Coefficient estimate") +
  scale_y_discrete(name="", labels=c("Time trend", "Concurrent\nnegotiation", "Post-\nnegotiation")) + 
  scale_x_continuous(breaks=breaks_pretty()) + 
  scale_color_manual("Negotiations", 
                     values=c("black", "gray30", "gray60")) +
  scale_shape_manual("Statistically\nsignificant (95%)", values=c(16, 15)) +
  facet_grid(cols=vars(fight), scales="free") + 
  theme(panel.spacing.x = unit(6, "mm"),
        legend.position = "bottom",
        legend.direction = "vertical")

ggsave("./figures/figure_a11.pdf", width=7, height=4.5)


##### APPENDIX 4.3: PLACEBO TESTS #####

# Make probabilities of picking wars
wn = unique(d$warName)
warp = NA
for (i in 1:length(wn)) {
  warp[i] = nrow(d[d$warName==wn[i],])/nrow(d)
}
wi = 60 # Temporal window
dsub = d |> dplyr::select(warName, date, neg, dayNum, all_of(fights))

cores = detectCores()
cl = makeCluster(cores[1]-4)
registerDoParallel(cl)


# Run the tests
afterEst = foreach (j = 1:1000, .combine='rbind') %dopar% {
  pl = list()
  
  set.seed(j)
  
  for (i in 1:150) {
    
    # Pick a random war
    war = sample(wn, 1, prob = warp)
    onewar = dsub[dsub$warName==war,]
    
    # Pick a random negotiation length
    nsim = sample(1:30, 1)
    
    # Pick a random day in the war; make negotiation start/end dates
    rday = sample(onewar$date, 1)
    nstart = rday
    nend = rday + nsim
    
    # Make sure dates are internal to the war
    nstart = max(head(onewar$date, 1) + 1, nstart)
    nend = min(tail(onewar$date, 1) - 1, nend)
    
    # Get the "wi" days before the start
    predays = onewar[which(onewar$date < nstart),]
    if (nrow(predays) > wi) { # Get only last part if too long
      predays = tail(predays, wi)
    }
    predays$afterNeg = 0
    
    # Get the "wi" days after the end
    postdays = onewar[which(onewar$date > nend),]
    if (nrow(postdays) > wi) { # Get only first part if too long
      postdays = head(postdays, wi)
    }
    
    postdays$afterNeg = 1
    
    # Put pieces together
    seg = rbind(predays, postdays)
    pl[[i]] = seg
  }
  tt = do.call("rbind", pl)
  
  # Now analyze the post-negotiation effect
  fit1 = lm(unlist(tt[,fights[1]]) ~ afterNeg + neg + dayNum + 
              factor(warName), data=tt)
  aos = coef(fit1)["afterNeg"]
  
  fit2 = lm(unlist(tt[,fights[2]]) ~ afterNeg + neg + dayNum + 
              factor(warName), data=tt)
  mom = coef(fit2)["afterNeg"]
  
  fit3 = lm(unlist(tt[,fights[3]]) ~ afterNeg + neg + dayNum + 
              factor(warName), data=tt)
  apd = coef(fit3)["afterNeg"]
  
  fit4 = lm(unlist(tt[,fights[4]]) ~ afterNeg + neg + dayNum + 
              factor(warName), data=tt)
  apddoe = coef(fit4)["afterNeg"]
  
  data.frame(aos, mom, apd, apddoe)
}
stopCluster(cl)

# Function to ensure parallel clustering is fully unregistered
unregister_dopar <- function() {
  env <- foreach:::.foreachGlobals
  rm(list=ls(name=env), pos=env)
}
unregister_dopar()

###### TABLE A26: 95% confidence intervals for placebo tests ######
round(summary(afterEst$aos), 3)
round(quantile(afterEst$aos, c(0.025, 0.975)), 3) # 95% CI: [-0.218, 0.392]

round(summary(afterEst$mom), 3)
round(quantile(afterEst$mom, c(0.025, 0.975)), 3) # 95% CI: [-0.378, 0.356]

round(summary(afterEst$apd), 3)
round(quantile(afterEst$apd, c(0.025, 0.975)), 3) # 95% CI: [-0.024, 0.008]


###### FIGURE A12: Bootstrapped estimates for the effect of the post-negotiation period on different battlefield measures ######

####### Figure A12(a): Recent imbalance #######
quickplot(afterEst$aos, geom="density") + 
  xlab("Coefficient estimate") + ylab("Density") + 
  geom_density(fill="gray80", color="gray80") + 
  geom_vline(xintercept=0, lty=2, linewidth=0.7) + 
  theme_bw() 
ggsave("./figures/figure_a12a.pdf", width=3.75, height=2.5, units="in")

####### Figure A12(b): Momentum #######
quickplot(afterEst$mom, geom="density") + 
  xlab("Coefficient estimate") + ylab("Density") + 
  geom_density(fill="gray80", color="gray80") + 
  geom_vline(xintercept=0, lty=2, linewidth=0.7) +
  theme_bw() 
ggsave("./figures/figure_a12b.pdf", width=3.75, height=2.5, units="in")

####### Figure A12(c): Inconsistency #######
quickplot(afterEst$apd, geom="density") + 
  xlab("Coefficient estimate") + ylab("Density") + 
  geom_density(fill="gray80", color="gray80") + 
  geom_vline(xintercept=0, lty=2, linewidth=0.7) + 
  theme_bw() 
ggsave("./figures/figure_a12c.pdf", width=3.75, height=2.5, units="in")


##### APPENDIX 4.4: ALTERNATIVE MEASURE OF CAPABILITIES #####

## Poisson regressions, using DOE scores instead of CINC scores ##
controlsB2 = "ceasefire + AOS60 + AOA + issues + contig + doeWinLose + polBinA + maj + nuk + 
  nBellig + keyAllies + anyDipEnemy + logBattles + nBattleTrend"

batFit1ad = glm(as.formula(paste("nBattle ~ neg + factor(warName)",
                                 controlsB2, sep="+")), data=d, family="poisson")
batFit1adr = coeftest(batFit1ad, vcov=vcovCL, cluster = ~ warName)

batFit1bd = glm(as.formula(paste("nBattle ~ negInt + negExt + factor(warName)",
                                 controlsB2, sep="+")), data=d, family = "poisson")
batFit1bdr = coeftest(batFit1bd, vcov=vcovCL, cluster = ~ warName)

car::linearHypothesis(batFit1bd, "negInt = negExt", singular.ok=T) # p << 0.01

# Change in number of battles
round(exp(batFit1ad$coefficients["neg"]), 3)    # 0.687
round(exp(batFit1bd$coefficients["negExt"]), 3) # 0.941
round(exp(batFit1bd$coefficients["negInt"]), 3) # 0.558


###### TABLE A27: Results of Poisson regressions, replacing CINC ratios with DOE scores ######

# Initial results before clustered standard errors
# Note that the number of observations per model only shows up here
stargazer(batFit1ad, batFit1bd,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum", "AOS60:logDayNum", "warName", "issues", "contig"),
          order = c("neg", "negInt", "negExt", 
                    "ceasefire", "AOS60", "AOA", "issues", "contig", "doeWinLose", "polBinA",
                    "maj", "nuk", "nBellig", "anyDipEnemy", "keyAllies", "logBattles"),
          covariate.labels = c("All negotiation", "Internal negotiation", "External negotiation",
                               "Ceasefire", "Recent imbalance", "Overall imbalance",
                               "DOE score",
                               "Democratic initiator", "Major power", "Nuclear", "Number of states",
                               "Opp. dip. representation",
                               "Major allies",
                               "Completed battles", "Active battle trend"),
          dep.var.labels.include=T,
          dep.var.labels = c("Active battles", "Active battles")
)

# Results after clustered standard errors
# Note that the number of observations per model does not show up here
stargazer(batFit1adr, batFit1bdr,
          column.sep.width = "0pt", no.space=T, align=T, df=F,
          omit = c("AOA:logDayNum", "logDayNum:AOA", "logBattles:logDayNum", "AOS60:logDayNum", "warName", "issues", "contig"),
          order = c("neg", "negInt", "negExt", 
                    "ceasefire", "AOS60", "AOA", "issues", "contig", "doeWinLose", "polBinA",
                    "maj", "nuk", "nBellig", "anyDipEnemy", "keyAllies", "logBattles"),
          covariate.labels = c("All negotiation", "Internal negotiation", "External negotiation",
                               "Ceasefire", "Recent imbalance", "Overall imbalance",
                               "DOE score",
                               "Democratic initiator", "Major power", "Nuclear", "Number of states",
                               "Opp. dip. representation",
                               "Major allies",
                               "Completed battles", "Active battle trend"),
          dep.var.labels.include=T,
          dep.var.labels = c("Active battles", "Active battles")
)


###### FIGURE A13: Coefficient plots for Poisson models of active battles with war fixed effects, replacing CINC ratios with DOE scores ######
battleCoefPlot(batFit1adr, batFit1bdr, "DOE")
file.rename("./figures/battleCoefPlotDOE.pdf",
            "./figures/figure_a13.pdf")

###### TABLE A28: Results of ordinary least squares regressions analyzing changes in battlefield measures, ######
# using all negotiations and a 60-day temporal window around failed negotiations 
stargazer(turn_all_aos_60, turn_all_mom_60, turn_all_apd_60, turn_all_apddoe_60,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"),
          dep.var.labels = c("Rec. imb.", "Momentum", "Inconsistency", "Incon. (DOE)"))

###### TABLE A29: Results of ordinary least squares regressions analyzing changes in battlefield measures, ######
# using the first negotiation in each war and a 60-day temporal window around failed negotiations 
stargazer(turn_early_aos_60, turn_early_mom_60, turn_early_apd_60, turn_early_apddoe_60,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))

###### TABLE A30: Results of ordinary least squares regressions analyzing changes in battlefield measures, ######
# using all negotiations after the first in each war and a 60-day temporal window around failed negotiations 
stargazer(turn_later_aos_60, turn_later_mom_60, turn_later_apd_60, turn_later_apddoe_60,
          column.sep.width = "-20pt", no.space=T, align=T, df=F,
          omit=c("warName"),
          covariate.labels = c("Post-negotiation", "Negotiation", "Time trend"))


###### FIGURE A14: Coefficient plots for OLS models of changes in battlefield outcomes after failed negotiations, using a 60-day temporal window ######
# Results for inconsistency calculated using DOE scores included.
revData |> 
  ggplot(aes(est, fct_rev(vars), group=fct_rev(model))) + 
  geom_vline(xintercept=0, linetype=2) + 
  geom_pointrange(aes(xmin=lower95, xmax=upper95, color=model, shape=signif), position=position_dodge(width=0.6)) + 
  theme_bw() + xlab("Coefficient estimate") +
  scale_y_discrete(name="", labels=c("Time trend", "Concurrent\nnegotiation", "Post-\nnegotiation")) + 
  scale_x_continuous(breaks=breaks_pretty()) + 
  scale_color_manual("Negotiations", 
                     values=c("black", "gray30", "gray60")) +
  scale_shape_manual("Statistically\nsignificant (95%)", values=c(16, 15)) +
  facet_grid(cols=vars(fight), scales="free") + 
  theme(panel.spacing.x = unit(6, "mm"),
        legend.position = "bottom",
        legend.direction = "vertical")
ggsave("./figures/figure_a14.pdf", width=7, height=4.5)